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This paper is a step towards a systematic theory of the transitivity (clustering) phenomenon 
in random networks. A static framework is used, with adjacency matrix playing the role of the 
dynamical variable. Hence, our model is a matrix model, where matrices are random, but their 
elements take values and 1 only. Confusion present in some papers where earlier attempts to 
incorporate transitivity in a similar framework have been made is hopefully dissipated. Inspired by 
more conventional matrix models, new analytic techniques to develop a static model with non-trivial 
clustering are introduced. Computer simulations complete the analytic discussion. 
PACS numbers: 02.50.Cw, 05.40.-a, 05.50.+q, 87.18.Sn 



I. INTRODUCTION 

Network model builders are currently adopting one of 
the two complementary approaches: Either a network is 
constructed step by step, by adding successive nodes and 
links. Or else, what is constructed is a static statisti- 
cal ensemble of networks. Each of these two approaches 
has its merits and shortcomings: Evolving network mod- 
els shed light on the growth dynamics, while static en- 
sembles are more appropriate for the study of structural 
traits. The classical model of Erdos and Renyi has 
been generalized so as to incorporate arbitrary degree 
distributions and even some correlations, but a serious 
shortcoming of the static models proposed so far is that 
they do not capture the common feature of most real 
networks: neighbors of a randomly chosen node are di- 
rectly linked to each other much more frequently than 
by chance, so that many short loops appear. The net- 
works tend to have locally a tree structure (see the re- 
views H H 3) And, as pointed out in ref. 0, "for 
general networks we currently have no idea how to in- 
corporate transitivity into random graph models" . In this 
paper we fill this gap, at least partially. The attention of 
the reader should be called to the very recent refs. 0-|f|. 
where the clustering problem is also addressed, but fol- 
lowing very different avenues. 

Graphs are a mathematical representation of networks. 
For definiteness we consider in this paper undirected 
graphs only. Let us denote by N the number of nodes in 
a graph and by M — {My}, i, j = 1, 2, TV the symmet- 
ric incidence matrix, with Ma = on the diagonal and 
Mij — 1 or depending on whether the nodes labeled 
i and j are connected or not. Whole information about 
a graph is encoded in its adjacency matrix. A general 
random graph model can be defined by introducing the 
partition function (sec, for example, ref Q): 



Z = J2 e S(JU) S(Tr(M 2 ) ~ 2L) 



(1) 
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where L is the number of links and S(M) is a function 
which we will call the action. The sum is over all possible 
adjacency matrices M. The simplest choice is S(M) = 
0. The corresponding graphs are those of the classical 
theory of Erdos and Renyi |]J: the value of the ratio 
L/N determines a variant of the model. 

Probably the simplest extension of the classical the- 
ory consists in setting S(M) — gTr(M 3 ), directly pro- 
portional to the number of triangles. This has been at- 
tempted already many years ago by Strauss @. His re- 
sults are summarized in the recent review 0: "There is 
however, one unfortunate pathology ... If, for example, 
we include a term in the Hamiltonian that is linear in 
the number of triangles in the graph, with an accompa- 
nying positive temperature favoring these triangles, then 
the model has the tendency to "condense" , forming re- 
gions of the graph that are essentially complete cliques - 
subset of vertices within which every possible link exists... 
Network in the real world however do not seem to have 
this sort of "clumpy" transitivity" . 

It appears to us that this negative conclusion, which 
faithfully reflects the content of ref. @, is not quite right. 
There is nothing wrong in Strauss's work. However, it is 
very incomplete and due to this incompleteness unvolun- 
tarily misleading. One of the aim of our paper is to give 
a fresh and comprehensive discussion of Strauss's model. 

The essence of Strauss's argument is as follows: assum- 
ing that the ratio L/N is kept constant, one can easily 
convince oneself that there exist pathological configura- 
tions for which Tr(M 3 ) oc iV 3 / 2 . The contribution of 
such a configuration to the partition function is explo- 
sive in the large N limit, since it cannot be tamed by the 
entropy factor falling roughly speaking like the inverse of 
the number of graphs, i.e. like exp (— const x NlogN). 
Thus, whatever small the coupling g is, the only stable 
states of the system are the pathological ones, provided 
the system is large enough. 

As we will show later on, the pathological crumpled 
states - the Strauss phase - are separated from a smooth 
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phase by a barrier that grows with increasing N. If the 
system is prepared in the smooth phase, it has a very 
tiny probability to roll out over the barrier to the Strauss 
phase. This probability tends rapidly to zero in the ther- 
modynamic limit. Strauss has missed this point, because 
the systems he simulated were too small to signal the rel- 
ative stability of the smooth phase. Now, for all practical 
purposes one can work in the smooth phase, ignoring the 
instability. This is what one does on many occasions 
in physics, in particular in the context of matrix models, 
where the instability also goes away when the matrix size 
tends to infinity. 

We have mentioned matrix models on purpose. The 
theory of random matrices is an important branch of sta- 
tistical physics, with applications ranging from nuclear 
physics to string theory. Some of the techniques devel- 
oped in this theory can be adapted to a study of the 
model defined by (JTJ. This is also a matrix model, albeit 
dealing with rather special matrices: in standard matrix 
models the matrix elements are continuous random vari- 
ables. 

The form of the partition function turns out to be 
very convenient for numerical simulations. In analytical 
calculations it will be convenient to use a slightly differ- 
ent formulation of the model, getting rid of the ^-function 
and allowing small fluctuations of the number of links L. 
The partition function Z will be, up to a factor, the av- 
erage of exp (S) in Erdos-Renyi theory. We first assume 
that a link is occupied with probability p. Hence, for 
given N and L the Erdos-Renyi weight is 

p L (l - p f{N-l)/2-L = (1 _ NJV(Ar-l)/a ( l _ y-L (2) 

p 

This primary weight is further multiplied by exp (S). In- 
serting l|2|l. integrating over L and neglecting an irrele- 
vant factor, we obtain the modified partition function 

Z = J2 ^P {~\ In {- l)7Y(Af 2 ) + S(M)) (3) 

In short, we have traded the 5-function for a Gaussian. 

In most of this paper we set S(M) — gTr(M 3 ), as 
in ref. Thus, formally and up to a rescaling of the 
dynamical variable the model looks like the much studied 
matrix model 

^matrix = / dM exp (- l -Tr(M 2 ) + gTr(M 3 )) (4) 

where one integrates over all possible symmetric N x N 
matrices. The difference is in the integration measure, 
which is discrete in J3J) and continuous 

dM = Y[ dM i:j (5) 
i<i 

in J!}. This difference is crucial, of course, but we would 
rather like to insist on the similarities between the two 



models. In any case, the example of the matrix model is 
for us a guide in our study. 

The plan of the paper is as follows: In Sect. II we 
develop exp (S) in powers of S and discuss the features 
of the perturbation series obtained by integrating term 
by term. In Sect. IIA we recall how the behavior of the 
perturbation series reflects the existence of an instabil- 
ity of the theory, by considering two examples. In Sect. 
IIB we introduce a helpful diagrammatic representation 
of the perturbative contributions to the partition func- 
tion. These diagrams are counted in Sect. IIC. It is 
argued that at finite N the perturbation series is patho- 
logical, indicating that nonperturbative phenomena are 
in action. However, keeping only the terms that are non- 
vanishing in the limit N — > oo one gets, like in the matrix 
model ijljl. a convergent series. This series is summed in 
Sect. IID. We obtain a simple analytic formula for the 
average number of triangles. We also show that the in- 
troduction of the interaction gTr(M 3 ) leaves the degree 
distribution unmodified. The nonperturbative dynamics 
is studied in detail in Sect. Ill, using the Monte Carlo 
technique of numerical simulation. In a range of model 
parameters we find a remarkable agreement between the 
data and the pertubative predictions, showing that the 
nonperturbative phenomena are negligible in this range. 
However, at large enough coupling strength the pertur- 
bation theory breaks down, as expected. The transition 
point has an interesting scaling with N. This enables us 
to define the model so as to get a non-trivial behavior of 
the clustering coefficient. In Sect. IV we discuss possible 
generalizations. This section contains also a summary of 
this work and a conclusion. 



II. PERTURBATION SERIES 

A. An analogy 

Before entering into the main discussion of our prob- 
lem let us consider an elementary example, to help those 
readers who are not conversant with field theoretic argu- 
ments. 

Consider the following integral 

I= y[£Jteef>(-*/*+9*/s-***) (6) 

where e is infinitesimal and has been introduced only in 
order to satisfy purists: I can be regarded as the partition 
function of a particle subject to the combined action of a 
potential and of a heat bath. Formally, the integrand re- 
sembles the summand in (|3J), except that the integration 
variable is here just a number. 

Consider a random walk in the potential given by the 
exponent in . Assume that in some initial moment the 
particle is located at x = 0. This is a metastable state. 
The particle eventually rolls over the barrier and reaches 
the deep minimum of the free energy at x f=a g/Ae. As 
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is well known, the lifetime r of the metastable state is 
given by the Arrhenius formula [9j: 

r » e^ 2 (7) 

The decay of the metastable configuration is a nonper- 
turbative phenomenon. The escape time has an essential 
singularity as a function of the coupling g. Of course, 
this nonperturbative phenomenon only occurs at nonzero 
temperature. When (3 = oo the particle stays forever in 
its initial position. Notice that the transition is more 
a cross-over than a genuine phase transition. It occurs 
when the exponent in Q is of order unity, but the value 
of g where the transition occurs may slightly depend on 
how the random walk is performed. 

In more complicated models the Arrhenius formula is 
not so readily derived. But nonperturbative dynamics 
shows up, if present, in the structure of the perturba- 
tion series in the coupling constant. Let us expand the 
exponential in @ with respect to terms other than the 
quadratic one: 

! r(3fc+l) 8. 9 2 fe 

1 ^ \ r(2fc + i) t W ( { )) { ) 

It is evident that the series coefficients grow factorially 
and that the series has zero radius of convergence. This 
is a characteristic signal. We will find a similar behavior, 
at finite N, in the model defined by ©. 

To conclude this section let us mention what happens 
when the partition function instead of being a Riemann 
integral, like in ©, is a matrix integral, like in 10}: A per- 
turbation expansion can again be defined and the terms 
in the expansion can be given a diagrammatic represen- 
tation. A clever manner of cataloging these diagrams 
has been devised by 't Hooft ^(J- One first rescales 
M -»■ \/NM,g -> g/\fN,e e/N. One then observes 
that these diagrams can be drawn on a two-dimensional 
surface. Such a surface is always a sphere with a num- 
ber of handles. Classes of diagrams are characterized by 
the number h of these handles and the contributions of 
all diagrams belonging to the same class have the same 
N dependence: N 2 ~ 2h . Summing over all h one gets a 
badly divergent series. However, in the limit N — > oo the 
spherical topology (h=0) dominates and the correspond- 
ing series has a non-zero radius of convergence. We will 
seek a similar behavior in our model. The hint is that 
one should carefully examine the N — > oo limit and that 
it may be wise to rescale the coupling constant in order 
to get a physically meaningful theory. 

B. Diagrammatic representation 

In this section we will introduce diagrams represent- 
ing terms in the perturbative expansion of To avoid 
misunderstanding let us stress from the outset that the 
diagrams introduced in this section are not to be identi- 
fied with the graphs belonging to the statistical ensemble 



we are working with. These diagrams are just a tool help- 
ing to catalogue contributions to the partition function. 
We start by setting S(M) = gTr(M 3 ). 
Let us expand in the factor e s : 

z= z °J2^(i Tr ( M3 w) ER 0) 

n 

Here Zq is the partition function in the Erdos-Renyi en- 
semble of random graphs and the subscript ER in (. . .) er 
indicates that the average is calculated in this ensemble. 
Since Zq does not depend on our dynamical coupling g 
it is for us an irrelevant normalization constant. 

The problem now is to calculate the averages appearing 
in the sum in |JS}. We use a method largely inspired by 
ref. adopting also some of their notations, like 

N& = Nl/(N -k)\ (10) 

which is the number of ways to choose k among N 
indices, the different permutations of the selected in- 
dices being considered as distinct. We have Tr(M 3 ) = 
J2abc MabMb c M ca , which is up to the factor 3! the num- 
ber T of triangles in the graph. We represent a matrix 
element M a b by a line segment. Indices a, b are then 
associated with the ends of the segment. The product 
M a bMi, c M ca is represented by a triangle. Notice, that 
M a t,Mi, c M ca is a random variable which can only take 
values or 1. Since the diagonal elements of M are by 
definition the probability that this product equals one 
is p 3 . There are N- configurations of indices which corre- 
spond to nonvanishing contributions to the sum. Hence 

(TrM 3 ) ER = p 3 N^ (11) 

The power of p is equal to the number of the sides of the 
triangle and the underlined power of N is the number 
of indices one is summing over. The contribution to the 
perturbation series is, of course, gp 3 N-. 

A complication arises as one goes to higher order of 
the perturbation theory. As one multiplies the traces 
one produces strings M a bMf, c . . . M e t where some pairs 
of indices, possibly interchanged like ab and ba, repeat 
themselves referring to the same element of the adjacency 
matrix. The corresponding probability factor is then p 
and not some power of p. One has to identify such pairs 
of indices. It is also necessary to identify the indepen- 
dent summation indices and to count their number. One 
has also to count in how many distinct combinations the 
independent indices can appear in the string. All this 
may seem a bit confusing and is best explained with an 
example. 

Consider the second order term in n — 2. One 
deals with the strings that have the following structure 

M ah M bc M ca M de M ef M fd (12) 

When one sums over indices there appear N- terms where 
the indices abedef are all distinct. Then, all the six index 
pairs are also distinct and the probability of individual 



(a) (b) (c) (d) 

FIG. 1: Diagrams representing 0(g 2 ) contributions to the partition function. 




(a) (b) 




(e) (f) 




(i) 0) 




(m) (n) 

FIG. 2: Diagrams representing 0(g [i ) 

strings is p 6 . We illustrate this situation by drawing two 
triangles (because there are two traces) that are non- 
overlapping, like in Fig. Id. The corresponding contri- 
bution to the perturbation series is 

fig Id - |-p 6 A^ 

Another possibility is that the indices def are identical 
to abc but appear in a different order. For each choice of 
abc there are 6 such possibilities and there are N- such 
choices. The probability of the string taking value 1 is p 3 . 
We illustrate this situation by drawing two overlapping 
triangles, as in Fig. la. Remember, that there are 6 
manners of putting one triangle on top of another. Hence 

fig la^6p 3 A^ 

Still another possibility is that two and only two in- 
dices are equal. They necessarily belong to two dis- 
tinct traces. Thus one of the indices def equals ei- 
ther a, or 6, or c: there are three possible choices. 
Let us take one of them, say a. Then there are three 
possible structures for the second trace in the product: 
M ae M ef M fa , M ea M af M fe , and M e fM fa M ae . Notice 
that e and / are dummy indices. We illustrate this situ- 
ation by drawing two triangles with one common vertex. 




(c) (d) 




(g) (h) 




(o) (P) 

contributions to the partition function. 

Remember that there are 3 manners of attaching a tri- 
angle to a vertex of another triangle, like in Fig. lc. On 
the whole there are 3x3 = 9 arrangements of the five 
independent indices and six distinct index pairs. Thus 

fig ic - lyVws 

Finally, two indices among def can be identi- 
cal to two indices among abc. There are three 
choices, let us take ab. Then there arc six pos- 
sible structures for the second trace in the prod- 
uct: M ab M be M ea , M ba M ae M ebl M ae M eb M ba , 
M ea M ab M bei M be M ea M ab and M eb M ba M ae . We illus- 
trate this situation by drawing two triangles with one 
common edge, like in Fig. lb. Remember that attach- 
ing a triangle to two vertices can be done is six manners. 
There are 3x6=18 arrangements of the 4 independent 
indices and five distinct pairs: 

fig lb->|^18p 8 JV± 

The game can be extended to higher order n, although 
the number of diagrams increases very fast. The third 
order diagrams are listed in Fig. 2. The general rule is 
that the power of p equals the number of triangle edges, 
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TABLE I: Contributions to the partition function correspond- 
ing to the diagrams of Fig. 2 ; the common factor g 3 /3! is 
omitted. 



(a): 36iVV 


(b): 324iVV 


(c): 216JVV 


(d): 1627VV 3 


(e): 648iV% 7 


(f): 108iV% 7 


(g): 324iV%> 8 


(h): 18JVV 


(i): 3247VV 


(j): 324iVV 


(k): 216iVV 


(1): 1627V V 


(m): 27iVV 


(n): 54iVV 


(o): 27iV%> 9 


(p): liVV 



follow the evolution of the number of triangle vertices 
during the recursive process, it is much more tedious to 
keep track of the number of triangle edges Thus, 
one can write a recursion equation for the sum 



W, 



(n) 



(14) 



that is for the total number of diagrams of order n, with 
k triangle vertices. This recursion relation, in essence, 
summarizes the rules listed in the preceding section: 



W, 



(n+l) 



= k(k- l){k-2)W } 



(„) 



3(k-l)(k 



3(k - 2)W<% + W, 



r(n) 
fc-3 



(15) 



The coefficients result from elementary combinatorics. 
The initial condition is wj^ — Sk3- The first two it- 
erations are listed below. One can check that the num- 
bers match those given in Figs. 1 and 2, provided the 
weights of diagrams with the same number of vertices 
are summed. 



the number of triangle vertices appears as the underlined 
power of N and, in the n* n order, there is a factor g n /n\. 
The determination of the number of independent index 
arrangements is rather tedious. The best way is to pro- 
ceed recursively. 

Using the diagrams one can actually forget about in- 
dices. It is sufficient to construct the diagrams of order n 
by adding one triangle to the diagrams of order n — 1 in 
all possible manners. One has to multiply the number of 
arrangements factor in the target diagram of order n — 1 
by the number of ways the new triangle can be attached 
to it. These numerical factors should be added when a 
given diagram of order n can be constructed from several 
diagram of order n — 1. We repeat again the rules: 

- free triangle: factor 1 

- triangle attached to one vertex: factor 3 

- triangle attached to a pair of vertices: factor 6 

- triangle attached to three vertices: factor 6 
The general structure of the perturbation series is 



^o = £f-r£^>> 



km l 



(13) 



The summation over k goes from 3 to 3n. The power m 
is always < k and m = k corresponds to diagrams where 
there is one or several groups of triangles lying one on 
top of another. Clearly, these are the only diagrams that 
survive in the limit N — > oo. 



C. Counting diagrams 

The quantity appearing in (jl"3"|) is the number of 

paths leading to a given diagram topology. In a sense it is 
the number of diagrams of that topology. We are able to 
determine it recursively, step by step, but we are unable 
to give a general formula for it. It is relatively easy to 



6VF 3 (1) = 6 



i (2) = 

wi 2) = 



18W. 



(i) 



18 



9PF 3 (1) = 9 



W, 



(2) 



ir. 



(i) 



w, 
u 
w, 
w, 



(3 



I'- 
ll 



6VF 



(2) 



= 24VF Z 



(2) 



(2) 



36 

f 18W; 

f 36VF 



(2) _ 



(2) 



540 

9W. 



(2) 



1242 



= 12fW ( 



(2) 



60W, 



(2) 



12 W, 



(2) 



IF 3 (2) = 882 



= 90W ( 
= 18VF 



(2) 
(2) 



15W, 



(2) 



w, 



(2) 



243 



W, 



(2) _ 



27 



(2) 



1 



One can easily see that the numbers given above agree 
with those presented in Table [I] For example, the mul- 
tiplicity of the diagrams (b) and (c) is, respectively, 324 

(3) 

and 216, which gives together W 4 = 540 as expected in 
the third order for the sum of diagrams occupying k = 4 
vertices. 

We can also estimate the number of diagrams 
in the large order of the expansion, that is for n — > oo. 
As argued, this number is expected to grow faster than 
a factorial, reflecting the non-perturbative transition to 
Strauss's phase. As we will show this is indeed the case. 
Define the polynomial function 



k 



(16) 



Multiplying both sides of eq. I|15fl by y k 3 and summing 
over k one obtains the following differential equation: 

W( n+1 \y) = y 3 (d y + l) 3 W (n Hy) (17) 



6 



The solution is 

W^(y) = [y 3 (d y + lf] n -l 
which can also be rewritten as 

W^(y) = e-y[y 3 d 3 y ] n 



(18) 



(19) 



Let us assume for a moment that W^ n '(y) grows with n 
less rapidly than (n!) K , with some fixed k, uniformly in y. 
It is then meaningful to introduce a generating function 



W{x,y) = Y,W {n) {y)x n /{n\T 



This function has a formal expansion 

>3 



z = er y yv 



where 



y ' ^ (n\) K 



(20) 



(21) 



(22) 



Developing e y in a power series we have 
Z -- 



y ^ W(xk(k-l)(k-2))y k 
"£ k\ ' (23) 



To check the convergence of this sum we need the asymp- 
totic behavior of a function W(z). Examples of such 
functions for integer n — 0,1,2 are well known, being 
a simple exponential e 2 , the Bessel function Iq(2^/z) or 
the generalized hypergeometric function (see |l3|) 0-^2(2) 
respectively. For arbitrary k one has 



W(z) 



C K 



Z 2 K 



(! + •••) 



(24) 



with some K-dependent constant C K . It is obvious that 
if k < 3 the series (|23|) is meaningful only when x = 0. 
For k > 3 it becomes convergent for arbitrary y and x. 
We conclude that W^ n ' grows faster than (n!) 3 ~ € , but 
slower than (n!) 3 for arbitrary small e. Such an explosive 
behavior of the number of perturbation theory diagrams 
is a signal that nonperturbative phenomena are present. 
Indeed, it means that the coefficients of high powers of 
g and TV -1 are increasing dramatically with the order 
of perturbation theory: what was assumed to be just a 
perturbation is in fact huge 0] ! 



D. Summation of leading diagrams 

We are interested in the limit N — > oo with pN = a = 
const. The structure of the perturbation expansion is 
given by eq. (jl3|l . As already mentioned, in general, the 
number of triangle edges (denoted m in I|1H|I) is larger 
or equal to the number of triangle vertices (denoted fc 
in 113(1 ). In the limit under consideration, only those 



diagrams contribute to the leading iV-independent term 
for which the number of triangle edges is equal to the 
number of triangle vertices. One can easily see that in 
these diagrams the triangles can overlap, but otherwise 
do not touch. In the expansion in g up to the third 
order, the following diagrams belong to this class: a single 
triangle in the first order, diagrams in Figs. la,d in the 
second, and those in Figs. 3 a,h,p in the third. Using the 
previously found results we have up to the third order: 

Z(G, 7 ) = 1+^7+^ (7 + 7 2 ) + ^ (7 + 3 7 2 + 7 3 )+- ■ • 

(25) 

where the convenient notation G = 6(7,7 = « 3 /6 has 
been introduced. In general, one can write this expansion 
as follows: 



G" 



G r - 



(n) fc 



^7) = l + £^(7) = l + £^£. fc '7 

n— 1 n— 1 k— 1 

(26) 

The coefficients z^ can be interpreted as the number 
of all diagrams which consist of n triangles located at 
k isolated positions, with possible multi-occupation of a 
position. Hence 



» 



£ 



(ni!) m i(n 2 !)" 12 . . . (n fc !) mfc mi!m 2 ! . . .m k \ 

(27) 

where n\ > n% > . . . n% and the sum is over all the par- 
titions P of n: nymx + niuii + . . . + nkniu = n. It turns 
out that the numbers z^ 1 can be calculated recursively: 

» , » 



4 n+1) =fe*r+4-'i fc = l,2,...,n + l (28) 



with the initial condition: z[^ — 1 and zi = for k 
outside of the closed interval [l,n]. The meaning of the 
equation is as follows. If one adds a new triangle to a 
configuration with n triangles, this triangle can be put 
at either of the k existing positions or at a new position. 
Hence, all configurations with n + 1 triangles located at 
k positions can be obtained from configurations with n 
triangles at k positions, by placing a new triangle at one 
of the k old positions, or from configurations with n tri- 
angles at [k — 1) positions by placing a new triangle at 
a new position. The first few terms resulting from this 
recursion relation are: 



z (2) h) = 


17 1 






z {2) {i) = 


17 1 - 


h l7 2 




^ (3) (7) = 


17 1 - 


h 3 7 2 4 


-I7 3 


^ (4) (7) = 


I7 1 - 


h 7 7 2 4 


- 6 7 3 + 7 4 


2 (5) (7) = 


17 1 - 


h 15 7 2 ■ 


+ 25 7 3 + IO7 4 + 7 5 


^ (6) (7) = 


17 1 - 


h 31 7 2 


+ 9O7 3 + 65 7 4 + 15 7 5 + 7 6 



The recursion relation can be converted into a partial 
differential equation for Z(G, 7). Multiplying both sides 



7 



of the equation by 7 and summing over k one finds 



<T> 



where 



i z (n+l) = jL Z (n) + z (n) 

7 dj 



fc=l 



(29) 



(30) 



Now, multiplying both sides by G n /n! and summing over 
n one obtains 



(31) 



where Z is given by (|26|) . An even simpler equation is 
satisfied by F = In Z: 



9 F ^ 7-! 1 



(32) 



One easily checks that the general solution is 

F(7,G) = /( 7 e G )-7 (33) 

where / is an arbitrary differentiable function. It results, 
however, from (|26|l that F{^, 0) = 0. Hence 



F(7,G)=lnZ( 7 ,G)= 7 (e G -l) 



(34) 



This result is not surprising for a practitioner of quantum 
field theory. Indeed, In Z(j, G) should equal the sum of 
contributions of connected diagrams. The only connected 
diagrams, in the large N limit, are those where triangles 
are all put on top of each other and according to our rules 
the diagram of n-th order yields just (g n /nl)p 3 Q n ~ 1 N- ~ 
jG n /nl. 

Notice, that the same result (|34|) is obtained assuming 
that in the Erdos-Renyi model the number of triangles 
has a Poisson distribution with average 7. Indeed, the 
average of exp (GT) is 



T=0 



E^e-V^exp^ 



1)] 



The average number of triangles is 



(T) 



d_ 

dG 



hxZ 



7e 



(35) 



(36) 



It is important to note that in the N — + 00 limit the 
average degree of a graph node becomes independent of 
G and is just equal to a, like in pure Erdos-Renyi theory. 
This can be easily seen in our formalism. Add to the 
action a source term r)Tr(PM), where P is the matrix 
Pij — Sn, so that Tr(PM) is the degree of the graph node 
with label 1. In our diagrammatics r/Tr(PM) produces 
a line instead of a triangle. But only one end of this line 
has a running index, the other end has index 1. The 
diagram of the lowest order in r\ is just this line and 




FIG. 3: The average number of triangles (T) versus the 
coupling constant G = 6g. The average degree is set to 
a — 4 and the simulation is performed for the number 
of nodes TV = 2 11 (squares), 2 13 (circles) and 2 15 (trian- 
gles) . The arrows indicate the position of the transition point 
G = Gout- The continuous line represents the analytic result 
<T) = (a 3 /6)exp(G). 



gives the contribution rjpN — rja. All corrections due 
to the interaction gTr(M 3 ) yield terms proportional to 
some inverse powers of N, because there is no way to 
put a triangle on a line. For example, the diagram of 
order i]g, where one has one triangle and one line on 
top of one of its edges gives 3rjgp 3 N- ~ 3r]ga 3 /N (we 
have N- and not N- because one of the triangle edges 
has the fixed label 1). In conclusion: the only connected 
diagram of order 0(rf) is independent of g, as is, to this 
order the free energy In Z. The average degree is just 
the derivative, at 77 = 0, of the free energy and equals a. 
One can extend this argument to higher order moments 
of the degree distribution. 

Using the results of this section we can propose a rough 
estimate of the expected region where nonperturbative 
physics sets in. With p = a/N and N large the summand 
in can be rewritten as 



N 

exp {In — (-L + G T)} 
a 



(37) 



We have rescaled the coupling by In ^ , so that this large 
factor multiplies now both terms in the action. We ex- 
pect that the perturbation series breaks down when the 
fluctuations of the two terms in the action become com- 
parable. The number of links is ~ N and we expect 
{(SL) 2 ) ~ N. In the large N limit the fluctuation of T 
is given by the second derivative of the free energy and 
equals ({ST} 2 ) = jexp(G) = jN G ". The two fluctua- 
tions become comparable when Go ~ 1. The numerical 
calculations confirm the logarithmic scaling of G, but as 
we will see in a moment the critical Go seems to lie below 
1. 



III. NUMERICAL SIMULATIONS 

We have shown in the preceding chapter that at finite 
N the perturbative series has a behavior which signals 



the presence of a nonperturbative phenomenon. The sim- 
ilarity of our problem with the example exhibited in the 
first subsection suggests an educated guess: there is a 
barrier separating the perturbative phase from a patho- 
logical but stable configuration; the nonperturbative phe- 
nomenon in question is the rolling of the system over 
the barrier towards this stable configuration. The bar- 
rier must become unpenetrable in the N — ► oo limit, be- 
cause in this limit the perturbation series becomes well 
behaved, actually summable, whatever is the coupling. 
We will confirm this guess with the help of numerical 
simulations. 





FIG. 4: As in Fig. EH (T) versus G, at N = 2 14 , but for three 
values of the average degree a — 2 (circles), 4 (squares) and 
8 (triangles). 



FIG. 5: The scaled clustering coefficient N{C) versus G at 
N = 2 14 and for a = 2 (circles), 4 (squares) and 8 (triangles). 
The line represents the expected behavior N{C) oc exp(G), 
where the proportionality coefficient is chosen so as to get the 
value expected in Erdos-Renyi model at G = 0. 




In constructing an algorithm manipulating adjacency 
matrices it is most important to take into account the 
sparse nature of these matrices. Only the positions of 
2L matrix elements carry a relevant information. This 
makes it possible to reduce the amount of computer mem- 
ory, needed to store an adjacency matrix, from 0(N ) in 
the naive coding to O(N) in the linear coding. In ef- 
fect, we simulate systems with the number of nodes of 
order 10 4 , i.e. three orders of magnitude larger then 
those simulated by Strauss |8(- In t he p resent work we 
use the algorithm introduced in refs. [lall3 by straight- 
forwardly upgrading it so as to include the term in the 
action proportional to the number of triangles. 

In the first numerical experiment we set a — 2L/N = 4 
and we measure the average number of triangles (T) for 
N = 2 11 , 2 13 and 2 15 . The coupling G is changed in small 
steps until the system makes a transition to Strauss's 
phase. The results are shown in Fig. |31 The continuous 
line corresponds to eq. (|36|l . It is remarkable that the 
points follow this line. The error bars are smaller than 
the symbol size. The transition points are indicated by 
an arrow. A closer examination of the data shows that 
near the transition the points start to deviate from the 
line and lie systematically above it. 



FIG. 6: The transition from the perturbative to Strauss's 
phase occurs at G — Gout- The figure shows how Gout de- 
pends on the system size N, for a = 2 (circles), 4 (squares) 
and 8 (triangles). Notice the logarithmic growth of the curves. 

In the next experiment we set N = 2 and measure 
(T) for a = 2, 4 and 8. The coupling G again varies up to 
the transition point. The result is shown in Fig. 0] The 
lines correspond to (T) = (a/6) exp (G). The agreement 
is remarkable. We have also measured the local clustering 
measure Cj as defined in ref. [la ]: 



where Tj is the number of triangles touching the vertex 
j, and Lj is the number of links emerging from it. We 
set Cj = when Lj is zero or one. A global clustering 
coefficient C is obtained by averaging over vertices. In 
Fig. 0we plot N(C) versus G. It is seen that N(C) = 
cr{a) exp (G), with a (a) = a(l — (1+a) exp (—a)) , which 
is the value of N(C) in Erdos-Renyi model. It is inter- 
esting that for very different values of a the transition 
occurs at roughly the same value of the clustering coeffi- 
cient. This is presumably not a numerical accident, but 
we have no explanation to offer. 
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FIG. 7: We compare the degree distribution calculated at two 
values of the coupling constant G. The number of nodes is set 
to N = 2048. At a — 2 the calculation is performed setting 
G = (triangles up) and 3.0 (circles). At a = 4 the calcula- 
tion corresponds to G — (squares) and G = 2.3 (triangles 
down). The continuous lines represent Poisson distributions 
with averages equal to 2 and 4, respectively. 



In Fig. El we show the variation of the transition point 
G = G ou t with the system size N . In our experiments 
the coupling G was always changed by SG — 0.1. Thus, 
in the figure we have associated an error 0.1 with the 
data points. Here a comment is in order: After having 
changed G we always made 1000 thermalization sweeps, 
then we carried out 20000 sweeps, performing measures 
every 10 sweeps. It is important to remember that the 
number of sweeps was always the same. Indeed, at finite 
N the system will sooner or later roll over the barrier, 
it is sufficient to wait long enough. The transition point 
G = G ou t is well defined when one decides to fix the 
waiting time. Actually, we are more interested by the 
scaling of G ou t with N than by its exact value. 

The curves in Fig. [13 are 

Go U t = 0.75 In AT - 2.4 for a = 2 (39) 
Gout = 0.70 In A" - 2.9 for a = 4 (40) 
Gout = 0.60 In A" - 2.7 for a = 8 (41) 

It is very interesting, although not really surprising (see 
Sect. IID), that G ou t scales like In AT. This means that 
setting G = Go In N one obtains a model with the clus- 
tering coefficient scaling non-trivially: C ~ AT 0-1 . 

Fig. illustrates the fact that the degree distribution 
is in the smooth phase insensitive to the value of the 
coupling G. We show the distribution at N = 2048 and 
a = 2 for G = and for a large value of G, i.e. G = 3.0, 
close to Gout ■ The distributions are almost identical and 
correspond to the Poisson distribution with average equal 
to a = 2 (the line) . This has been repeated for A" = 2048 
and a = 4, where we measured at G — and G = 2.3. 



IV. DISCUSSION, SUMMARY AND 
CONCLUSION 

A. Possible generalizations 

Up to now we assumed that the interaction has a 
simple form S(M) = gTr(M 3 ). The question which 
immediately comes to mind is: what happens to the 
network transitivity when the action is more compli- 
cated? Assume that S(M) has the polynomial form: 
S(M) = E„> 3 g n Tr{M n ), with 53 = g. Our diagram- 
matic rules can easily be extended to include this case. 
Tr(M n ) is represented by a polygon with n sides. How- 
ever the polygon can be folded, the same line segment 
being covered several times. In particular, when n > 4, 
a polygon can be folded so that some of its edges form a 
triangle. 

The calculation of leading diagrams given in Sect. IID 
can be generalized to interactions involving odd powers 
of M. One can limit oneself to connected diagrams, those 
contributing to the free energy. In the N — > oo limit the 
contribution of a diagram is proportional to A^ to a power 
equal to # vertices - # edges (because p = a/N). We are 
interested in diagrams which overlap with a triangle. In 
this case # vertices - # edges = 0. In all other cases this 
quantity is negative and the diagram does not contribute 
in the limit. It is easy to see that diagrams surviving in 
the A^ — > oo limit are those where triangles and folded 
polygons are put on top of each other. 

The situation is more complicated for even powers of 
M. The leading diagrams look like branched polymers, 
with # vertices = # edges + 1, and their contribution 
diverges like N. In order to avoid an unwanted renor- 
malization of the quadratic term in the action one has to 
subtract from Tr(M 2k ) a counterterm ~ Tr(M 2 ) with 
an appropriate coefficient in front. Then the calculation 
is like for the odd power case. 

We have not pushed this calculation very far. As far as 
we can see one expects a certain degree of universality: 
the higher powers of M should not change the qualitative 
picture very much, although they may be important for 
phenomenology, to fit the data. A comprehensive study 
of these more general interactions is certainly worth being 
done. This is, however, beyond the scope of the present 
paper. 

It is not quite clear what is the best way of extend- 
ing the theory of this paper so as to obtain an arbitrary 
degree distribution. The field theoretical methods ex- 
tensively used in this paper usually fail when the action 
becomes nonanalytic. The simplest, although perhaps 
not the most elegant, extension consists in using instead 
of the Erdos-Renyi model, a general model with uncorre- 
cted nodes |17| as the zeroth order approximation. Pre- 
liminary numerical results look encouraging, although it 
is clear that much has to be done in order to get a fully 
satisfactory phenomenology. We hope to return to this 
problem elsewhere. 



10 



B. Summary 

Let us now summarize what has been achieved in this 
work: In most of this paper we have discussed a model 
of random graphs where the classical Erdos-Renyi theory 
is generalized by the introduction into the action of an 
interaction term ^Tr(M 3 ) = GT, M being the N x N 
adjacency matrix and T the number of triangles, respec- 
tively. This model is our guinea pig. It is a matrix model, 
but of a special kind, because the dynamical variable is 
a random matrix whose elements equal either or 1. 

Inspired by the analogy with more conventional ma- 
trix models we develop a diagrammatic technique, en- 
abling us to calculate the perturbation series analytically. 
We count the diagrams and show that, at finite N, their 
number grows so rapidly that the perturbation series be- 
comes pathological. This also happens in conventional 
matrix models and, as is well known, indicates the pres- 
ence of a nonperturbative phenomenon. The nature of 
this phenomenon is identified through numerical simu- 
lations. There is a "potential barrier" and the system 
can roll over it and fall into a pathological phase, where 
all triangles form a unique clan. This phase was first 
discovered long ago by Strauss 0, who did not notice, 
however, that it is separated from a smooth phase by a 
barrier which becomes impenetrable at large N. We pro- 
pose to consider this smooth phase as the physical one. 

We show, that for large enough N and in a range of 
values of the coupling constant G the smooth phase can 
be considered, for all practical purposes, as stable. In this 
range of G it is meaningful to neglect the nonperturbative 
physics and to limit oneself to leading diagrams (those 
obtained setting N = oo). We are able to sum all these 
diagrams up, obtaining simple analytic expressions for 
the free energy and for the average number of triangles. 
We also show analytically that in this regime the degree 
distribution is insensitive to the value of the coupling G. 

A heuristic argument, confirmed by numerical simula- 
tions, indicates that the transition point G = G out , where 
the system jumps to the Strauss's phase, scales with N 



like In AT. Hence, the physical coupling is not so much G 
but rather Go defined by the equation G = Go In A. Our 
simulations indicate that at the transition point Go = 0.6 
to 0.75, depending on the average degree, but this result 
should be taken with a grain of salt. Anyhow, the cluster- 
ing coefficient scales non-trivially, like G ~ N 00 ^ 1 and 
is larger by one to two orders of magnitude than in the 
unperturbed model. 

It appears that the analytic treatment can be extended 
to more complicated, but polynomial actions. In the 
present state of affairs the extension of our approach to 
more realistic, for example scale-free models can only be 
done numerically. 

C. Conclusion 

Clustering is a rather striking trait of many observed 
networks. The local tree-like structure characterizing 
most static models is clearly non-realistic. We have ar- 
gued elsewhere that static models are an important in- 
gredient of network theory. Thus, we believe that it is 
important to be able to construct static models with non- 
trivial clustering. There was some confusion concerning 
the feasibility of such an enterprise. We hope having dis- 
sipated it. For the sake of clarity we have focused our 
attention on a model where much can be done analyti- 
cally. It is a specific matrix model, where matrices are 
random, but their elements take values and 1 only. In 
the zero-th order approximation it is equivalent to the 
classical Erdos-Renyi model of graphs. Non-trivial clus- 
tering is generated by an appropriate interaction. A com- 
prehensive phenomenologically oriented study is beyond 
the scope of this paper and remains to be carried out. 
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